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Abstract 

We give a short and elementary introduction to Lie group methods. A family 
of symplectic integrators on cotangent bundles of Lie groups is presented and the 
notion of discrete gradient methods is generalized to Lie groups. Finally a selection of 
applications of Lie group integrators are discussed. 

1 Introduction 

The significance of the geometry of differential equations was well understood already in 
the nineteenth century, and in the last few decades such aspects have played an increasing 
role in numerical methods for differential equations. Nowadays, there is a rich selection 
of integrators which preserve properties like symplecticity, reversibility, phase volume and 
first integrals, either exactly or approximately over long times [28]. Differential equations 
are inherently connected to Lie groups, and in fact one often sees applications in which the 
phase space is a Lie group or a manifold with a Lie group action. In the early nineties, 
two important papers appeared which used the Lie group structure directly as a building 
block in the numerical methods. Crouch and Grossman |20j suggested to advance the 
numerical solution by computing flows of vector fields in some Lie algebra. Lewis and 
Simo [12] wrote an influential paper on the Lie group based integrators for Hamiltonian 
problems, considering the preservation of symplecticity, momentum and energy. These 
ideas were developed in a systematic way throughout the nineties by several authors. In 
a series of three papers, Munthe-Kaas |501 [5Tj [52] presented what are now known as the 
Runge-Kutta-Munthe-Kaas methods. By the turn of the millennium, a survey paper [33] 
summarized most of what was known by then about Lie group integrators. 

The purpose of the present paper is three- fold. First, in section [2] we would like to give 
an elementary, geometric introduction to the ideas behind Lie group integrators. Secondly, 
we present some new material in sections [3] and [4] Symplectic Lie group integrators have 
been known for some time, derived by Marsden and coauthors [35] by means of variational 
principles. Here we consider a group structure on the cotangent bundle of a Lie group 
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and derive symplectic Lie group integrators using the model for vector fields on manifolds 
denned by Munthe-Kaas in [52] . In section [4] we extend the notion of discrete gradient 
methods as proposed by Gonzalez [27] to Lie groups, and thereby we obtain a general 
method for preserving first integrals in differential equations on Lie groups. Our third 
objective is to present a selection of applications of Lie group integrators. There are many 
such examples to choose from, and we give here only a few teasers. 

We would also like to briefly mention some of the issues we are not pursuing in this 
article. One is the important family of Lie group integrators for problems of linear type, 
including methods based on the Magnus and Fer expansions. An excellent review of the 
history, theory and applications of such integrators can be found in [2j. We will also skip 
all discussions of order analysis of Lie group integrators. This is a large area by itself which 
involves technical tools and mathematical theory which we do not wish to include in this 
relatively elementary exposition. There have been several new developments in this area 
recently, in particular by Lundervold and Munthe-Kaas, see e.g. [Hj. 

2 Lie group integrators 

The simplest consistent method for solving ordinary differential equations is the Euler 
method. For an initial value problem of the form 

y = F(y), y(0) = yo, 

one takes a small time increment h, and approximates y(h) by the simple formula 

yi = yo + hF(y ), 

advancing along the straight line coinciding with the tangent at yo- Another way of thinking 
about the Euler method is to consider the constant vector field F y o(y) := -F(yo) obtained 
by parallel translating the vector F(yo) to all points of phase space. A step of the Euler 
method is nothing else than computing the exact /i-flow of this simple vector field starting 
at yo- in Lie group integrators, the same principle is used, but allowing for more advanced 
vector fields than the constant ones. A Lie group generalization of the Euler method is 
called the Lie-Euler method, and we shall illustrate its use through an example [20]. 

Example 1, the Duffing equation. Consider the system in R 2 



a model used to describe the buckling of an elastic beam. Locally, near a point (xo, yo) we 
could use the approximate system 



x = y, 



a > 0,6 > 



(1) 



y = —ax 



x = y, 

y = -{a + bxl 



x(0) = x , 
2/(0) = yo, 



(2) 
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Figure 1: (R d , +)-frozen vector field (left) and SX(2)-frozen vector field (right) for the 
Duffing equation 



which has the exact solution 

x(t) = xq cosut + — sinwt, y(t) = yo cosuit — ojxq sinwt, 
Alternatively, we may consider the local problem 



\Ja + bx%. 



x 

y 



= —ax — bxy, 



(3) 



having exact solution 
x(t) 

m 



x cosat+fsmat + bxf, ^g" 1 
yo cos at — axo sin at 



jj x 3 sin at 
a 



a 



(4) 



In each of the two cases, one may take x\ = x(h), y\ = y(h) as the numerical approximation 
at time t = h. The same procedure is repeated in subsequent steps. A common framework 
for discussing these two cases is provided by the use of frames, i.e. a set of of vector fields 
which at each point is spanning the tangent space. In the first case, the numerical method 
applies the frame 



X 



V 




=: y dx, Y 





x 



-: x dy. 



(5) 



Taking the negative Jacobi-Lie bracket between X and Y yields the third element of the 
standard basis for the Lie algebra sl(2), i.e. 



H = — [X, Y] = x dx — y dy, 



(6) 
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so that the frame may be augmented to consist of {X, Y, H}. In the second case, the 
vector fields E\ = y dx — ax dy and E2 = dy can be used as a frame, but again we choose 
to augment these two fields with the commutator E3 = — \E\,E2[ = dx to obtain the Lie 
algebra of the special Euclidean group SE{2) consisting of translations and rotations in 
the plane. 

In general, a way to think about Lie group integrators is that we have a manifold M 
where there is such a frame available; {Ei, . . . , Ed} such that at any point p G M one has 

span{E 1 {p),...,E d {p)} = T p M 

and frames with this property are said to be locally transitive. This may be a linear space 
or in many cases even a Lie algebra of vector fields. In the example with Duffing's equation, 
the set {X,Y,H} is locally transitive on R 2 \{0} and {E±, E2, £3} is locally transitive on 
R 2 . 

Given an arbitrary vector field F on M, then at any point p £ M there exists a vector 
field F p in the span of the frame vector fields such that F p (p) = F(p). An explicit way of 
writing this is by using a set of basis vector fields £1,...,^ forg, such that any smooth 
vector field F has a representation 

d 

F(y) = Y J fk{y)E k {y), (7) 

k=l 

for some functions ff. : M — > R. The vector fields F p £ q, called vector fields with frozen 
coefficients by Crouch and Grossman |20| . are then obtained as 

d 

F P (y) = Y J fk(p)E k (y)- 

k=l 

In the example with the Duffing equation we took E\ = X, E2 = Y, fi(x,y) = 1 and 
f2{x, y) = —{a + bx 2 ). The Lie-Euler method reads in general 

y n+1 = exjp(hF yn )y n , 

where exp denotes the flow of a vector field. 

A more interesting example, also found in [20] is obtained by choosing M = S 2 , the 
2-sphere. A suitable way to induce movements of the sphere is that of rotations, that is, by 
introducing the Lie group SO(3) consisting of orthogonal matrices with unit determinant. 
The corresponding Lie algebra so (3) of vector fields are spanned by 

Ei(x,y, z) = —zdy + ydz, E2(x,y, z) = zdx — xdz, E^(x,y, z) = —ydx + xdy. 

We note that xE\(x, y, z) + yE2(x, y, z) + zE^(x, y, z) = 0, showing that the functions fk 
in (TF| are not unique. A famous example of a system whose solution evolves on S 2 is the 
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free rigid body Euler equations 



,1 1, 11 . 1 1 

T - T )yz, y=(—-—)xz, z={—-—)xy, 



(8) 



x 
h 





-z 

y 



y_ 
h 



z 


-X 



z 

h 



-y 

X 





where x,y,z are the coordinates of the angular momentum relative to the body, and 
Ii,l2, 13 are the principal moments of inertia. A choice of representation ([7]) is obtained 
with 

fl(x,y,z) = ~, f 2 (x,y,z) = -—, f 3 (x, y, z) = -—, 

+1 ±2 -*3 

so that the ODE vector field can be expressed in the form 
F(x,y,z) 

We compute the vector field with coefficients frozen at po = (xo, yo,zo), 

x 

y 

z 

The /i-flow of this vector field is the solution of a linear system of ODEs and can be 
expressed in terms of the matrix exponential expm(/ii ? p o). The Lie-Euler method can be 
expressed as follows: 








I 3 


yo 
h 


F po {x,y,z) = 


h 





.t'o 
h 




yo 

. h 


x 

h 






Given p = (x ,y ,z ) 
for n = 0, 1, . . . 

p n+ i = expm(/iF p J -p n 
end for 

Notice that the matrix to be exponentiated belongs to the matrix group so (3) of real 
skew-symmetric matrices. The celebrated Rodrigues' formula 

... sin a , 1 — cos a <2 2 ,, ,,,2 1,, ,,,2 , , nS 

expm(A) = I-\ A-\ ^ A , a = \\A \\i = \ \\A i, A £ so(3), 

a a z 

provides an inexpensive way to compute this. 

Whereas the notion of frames was used by Crouch and Grossman in their pioneering 
work [20], a different type of notation was used in a series of papers by Munthe-Kaas 
|50l [5Tj I52j . see also |44j for a more modern treatment. Let G be a finite dimensional Lie 
group acting transitively on a manifold M. This means that for any two points mi, rri2 G M 
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there exists a group element g G G such that 771,2 = 9 ■ mi- We denote the Lie algebra of G 
by g. For any element £ G g there exists a vector field on M 



d 

d7 



exp(t£) • 777 =: A*(£)(m) 



0) 



t=o 



Munthe-Kaas introduced a generic representation of a vector field F S Af(M) by a map 
f-.M—^Q such that 

F(m) = A*(/(m))(m) (10) 

The corresponding frame is obtained as Ei = A*(ej) where ei, . . . , is some basis for g 
and one chooses the functions /j: M — > R such that /(m) = £)i=i fi( m ) e i- The map A* 
is an anti-homomorphism of the Lie algebra g into the Lie algebra of vector fields X{M) 
under the Jacobi-Lie bracket, meaning that 



^*([X m ,Y m ] s ) — — [\*(X m ), A*(y^ 7l )] 



JL 



This separation of the Lie algebra g from the manifold M allows for more flexibility in 
the way we represent the frame vector fields. For instance, in the example with Duffing's 
equation and the use of sl(2), we could have used the matrix Lie algebra with basis elements 



An 



1 




Y 

1 rii 





1 



1 






-1 



rather than the basis of vector fields ([5]), ([6]). The group action by g G SL(2) on a point 
m £ R 2 would then be simply g • 777, matrix- vector multiplication, and the exp in Q would 
be the matrix exponential. The map f(x, y) would in our case be 





-(a + bx 2 ) 



y 





but note that since the dimension of the manifold is just two whereas the dimension of 
sl(2) is three, there is freedom in the choice of /. In the example we chose not to use the 
third basis element H. 



2.1 Generalizing Runge— Kutta methods 

In order to construct general schemes, as for instance a Lie group version of the Runge- 
Kutta methods, one needs to introduce intermediate stage values. This can be achieved in 
a number of different ways. They all have in common that when the methods are applied 
in the Euclidean space where the Lie group is (R m , +), they reduce to conventional Runge- 
Kutta schemes. Let us begin by studying the simple second order Heun method, sometimes 
called the improved Euler method. 

h = F(y n ), k 2 = F(y n + hki), y n +i = Vn + \ h(fa + k 2 ) 
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Geometrically, we may think of k\ and k% as constant vector fields, coinciding with the 
exact ODE F(y) at the points y n and y n + hk\ respectively. The update y n +i can be 
interpreted at least in three different ways, 

ex P(!|(^i + k 2 ))y n , exp(^fci) o exp(^k 2 )y n , exp(^/c 2 ) ° exp(^ki)y n (If) 

The first is an example of a Runge-Kutta-Munthe-Kaas method and the second is an 
example of a Crouch-Grossman method. All three fit into the framework of commutator- 
free Lie group methods. All three suggestions above are generalizations that will reduce to 
Heun's method in (R m +), in principle we could extend the idea to Runge-Kutta methods 
with several stages 

s s 

y n +\ =y n + h'^2 biF(Yi), Yi = y n + ft a{jF(Yj), i = l,...,s, 
i=i j=i 

by for instance interpreting the summed expressions as vector fields with frozen coefficients 
whose flows we apply to the point y n S M. But it is unfortunately not true that one in 
this way will retain the order of the Runge-Kutta method when applied to cases where the 
acting group is nonabelian. 

Let us first describe methods as proposed by Munthe-Kaas [52], where one may think 
of the method simply as a change of variable. As before, we assume that the action of G on 
M is locally transitive. Since the exponential mapping is a local diffeomorphism in some 
open set containg € 0, it is possible to represent any smooth curve y{t) on M in some 
neighborhood of a point p G M by means of a curve a(t) through the origin of g as follows 

y(t) = exp(a(t)) ■ p, a(0) = 0, (12) 

though a(t) is not necessarily unique. We may differentiate this curve with respect to t to 
obtain, 

y(t) = \*(dexp ait) d-(t))(y(t)) = F(y(t)) = A*(/(exp(a(t)) • p)(y(t)) (13) 

The details are given in [52j and the map dexp CT : g — > q was derived by Hausdorff in [30] 
as an infinite series of commutators 

1 r^..i , 1 r_r_„.n , _ V" 1 „jfc„. _ ex P( z ) ~ 1 



dexp CT (w) = v + -[a,v} + -[a, [a, v]} + --- = ^ ( fc + 1 ); ad ^ 



v (14) 

with the usual definition of &d u (v) as the commutator [u,v]. The map A* does not have to 
be injective, but a sufficient condition for (13) to hold is that 



a = dexp (T 1 (/(exp(cr) ■ p)) 

This is a differential equation for a(t) on a linear space, and one may choose any conven- 
tional integrator for solving it. The map dexp" 1 : g — > g is the inverse of dex.p a and can 
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also be interpreted as the right trivialized derivative of the logarithm. From (14) we find 
that one can write dexp~ 1 (u) as 



dexp^ 1 ^) 



exp(z) — 1 



v = v--[a,v} + — [a,[a, «]] + ••• (15) 



The coefficients appearing in this expansion are scaled Bernoulli numbers and i?2fc+i = 
for all /c > 1. One step of the resulting Runge-Kutta-Munthe-Kaas method is then 
expressed in terms of evaluations of the map / as follows 

s 

yi = expih^biki) ■ y 
i=i 



This is not so surprising seen from the perspective of the first alternative in (11), the main 
difference is that the stages ki corresponding to the frozen vector fields \*(ki) need to be 
"corrected" by the dexp -1 map. Including this map in computational algorithms may seem 



awkward, however, fortunately truncated versions of (15) may frequently be used. In fact, 
by applying some clever tricks involving graded free Lie algebras, one can in many cases 
replace dexp -1 with a low order Lie polynomial while retaining the convergence order of 
the original Runge-Kutta method. Details of this can be found in |53| 15]. There are also 
some important cases of Lie algebras for which dexp" 1 can be computed exactly in terms of 
elementary functions, among those is so (3) reported in [15]. Notice that the representation 



(12) does not depend on the use of the exponential map from g to G. In principle, one can 
replace this map with any local diffeomorphism ip, usually one scales ip such that <^(0) = e 
and TQip = Id g . Examples of such maps are the Cayley transformation |23| which can be 
used for matrix Lie groups of the type Gp = {X 6 j^ dxd : X T PX = P} for a nonsingular 
d x d- matrix P. These include the orthogonal group 0(n) = Gj and the linear symplectic 
group SP(n) = Gj where J the skew-symmetric matrix of the standard symplectic form. 
Another possibility is to replace the exponential map by canonical coordinates of the second 
kind [53]. 

We present here the well-known Runge-Kutta-Munthe-Kaas method based on the pop- 



S 



ular fourth order method of Kutta [36], having Butcher tableau 
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1 1 




6 


3 


3 6 



(16) 



In the Lie group method, the dexp 1 map has been replaced by the optimal Lie polynomials. 
h = hf(y ), 
k 2 = /i/(exp(|fci) • y ), 
h = fo/(exp(±& 2 - l[h,k 2 ]) ■ y ), 
k A = /i/(exp(fc 3 ) • y ), 

yi = exp(|(A?i + 2k 2 + 2k 3 + k± - § [fci , fc 4 ] ) ) • y Q . 

An important advantage of the Runge-Kutta-Munthe-Kaas schemes is that it is easy to 
preserve the convergence order when extending them to Lie group integrators. This is 
not the case with for instance the schemes of Crouch and Grossman [201 l5Tj , where it is 
necessary to develop order conditions for the nonabelian case. This is also true for the 
commutator- free methods developed by Celledoni et al. [13] . In fact these methods include 
those of Crouch and Grossman. The idea here is to allow compositions of exponentials 
or flows instead of commutator corrections. With stages k\, . . . ,k s in the Lie algebra, one 
includes expressions of the form 

exp(J^ ft!*) ■ ■ ■ exp(^ [3 l 2 ki) exp(^ (3\ki)y Q 



both in the definition of the stages and the update itself. In some cases it is also possible 
to reuse flow calculations from one stage to another, and thereby lower the computational 
cost of the scheme. An extension of (16) can be obtained as follows, setting k{ = hf(Yi) 
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for all i, 

Yi = yo, 

Y 2 = exp(|/ci)y , 

Y 3 = exp( \k 2 )y 

Y i = eMKh-lh))Y 2 , 

yi = exp(^(3A;i + 2k 2 + 2k 3 - fc 4 )) y , 

yi = exp(^(-fci + 2k 2 + 2k 3 + 3fc 4 )) yi. 

Note in particular in this example how the expression for I4 involves Y 2 and thereby one 
exponential calculation has been saved. 

2.2 A plenitude of group actions 

We saw in the first examples with Duffing's equation that the manifold M, the group G 
and even the way G acts on M can be chosen in different ways. It is not obvious which 
action is the best or suits the purpose in the problem at hand. Most examples we know 
from the literature are using matrix Lie groups G C GL(n), but the choice of group action 
depends on the problem at hand and the objectives of the simulation. We give here several 
examples of situations where Lie group integrators can be used. 

G acting on G. In the case, M = G it is natural to use either left or right multiplication 
as the action 

L g (m) = g • m or R g (m) = m ■ g, g,m € G. 



The correspondence between the vector field F 6 X(M) and the map (10) is then just the 
tangent map of left or right multiplication 

F(g) = T e L g (f(g)) or F(g) = T e R g (f(g)), g G G. 

When working with matrices, this simply amounts to setting F{g) = g ■ f(g) or F{g) = 
f(g) ■ g. Note that f(g) is related to f(g) through the adjoint representation of G, Ad : 
G -> Aut(fl) 

f(g) = Ad g f(g), Ad g = T e L g o T e R~\ 
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The coadjoint action and Lie Poisson systems. Lie group integrators for this in- 
teresting case were studied by Eng0 and Faltinsen [25] . Suppose G is a Lie group and the 
manifold under consideration is the dual space q* of its Lie algebra q. The coadjoint action 
by G on g* is denoted Ad* defined for any g G G as follows: 

(Ad>,0 = </i,Ad ff £), V£e . 

It is well known that mechanical systems formulated on the cotangent bundle T*G with a 
left or right invariant Hamiltonian can be reduced to a system on q* given as 

/i = ±ad%H H 

the negative sign is used in case of right invariance. The solution to this system preserves 
coadjoint orbits, which makes it natural to suggest the group action 

9 ■ M = A d*_i^, 

so that the resulting Lie group integrator also respects this invariant. For Euler's equations 
for the free rigid body, the Hamiltonian is left invariant and the coadjoint orbits are spheres 
in g* ^ R 3 . 

Homogeneous spaces and the Stiefel and Grassmann manifolds. The situation 
when G acts on itself by left of right multiplication is a special case of a homogeneous 
space [55] , where the assumption is only that G acts transitively and continuously on some 
manifold M. Homogenous spaces are isomorphic to the quotient G/G x where G x is the 
isotropy group for the action at an arbitrarily chosen point x £ M 

G x = {y £ G : y ■ x = x}. 

All isotropy groups G x , x £ M are conjugate and therefore the isomorphism type of 
G/G x is independent of x. A much encountered example is the hypersphere M = S^ 1 
corresponding to the left action by G = SO(d), the Lie group of orthogonal d x d matrices 
with unit determinant. One has S^ 1 = SO{d) / SO(d — 1). We have in fact already 
discussed the example of the free rigid body ^ where M = S 2 . 

The Stiefel manifold St(d, k) can be represented by the set of d x £>matrices with 
orthonormal columns. An action on this set is obtained by left multiplication by G = 
SO(d). Lie group integrators for Stiefel manifolds are extensively studied in the literature, 
see e.g. [HI [35] and some applications involving Stiefel manifolds are discussed in Section[6j 
An important subclass of the homogeneous spaces is the symmetric spaces, also obtained 
through a transitive action by a Lie group G, where M = G/G x , but here one requires 
in addition that the isotropy subgroup is an open subgroup of the fixed point set of an 
involution of G [54J. A prominent example of a symmetric space in applications is the 
Grassmann manifold, obtained as SO(d)/(SO(k) x SO(d — k)). 
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Isospectral flows. In isospectral integration one considers dynamical systems evolving 
on the manifold of d x d-matrices sharing the same Jordan form. Considering the case of 
symmetric matrices, one can use the transitive group action by SO(d) given as 

T 

g ■ m = g mg. 

This action is transitive, since any symmetric matrix can be diagonalized by an appropri- 
ately chosen orthogonal matrix. If the eigenvalues are distinct, then the isotropy group is 
discrete and consists of all matrices in SO(d) which are diagonal. 

Lie group integrators for isospectral flows have been extensively studied, see for example 
[U [7] . See also [9] for an application to the KdV equation. 



Tangent and cotangent bundles. For mechanical systems the natural phase space will 
often be the tangent bundle TM as in the Lagrangian framework or the cotangent bundle 
T*M in the Hamiltonian framework. The seminal paper by Lewis and Simo [12] discusses 
several Lie group integrators for mechanical systems on cotangent bundles, deriving meth- 
ods which are symplectic, energy and momentum preserving. Eng0 [21] suggested a way to 
generalize the Runge-Kutta-Munthe-Kaas methods into a partitioned version when M is 
a Lie group. Marsden and collaborators have developed the theory of Lie group integrators 
from the variational viewpoint over the last two decades, for an overview, see |45| . For 
more recent work pertaining to Lie groups in particular, see |38|. El [61]. In Section [3] we 
present what we believe to be the first symplectic partitioned Lie group integrators on 
T*G phrased in the framework we have discussed here. Considering trivialized cotangent 
bundles over Lie groups is particularly attractive since there is a natural way to extend 
action by left multiplication from G to G x q* via (19). 



The affine group and its use in semi-linear PDE methods. Lie group integrators 
can also be used for approximating the solution to partial differential equations, the most 
obvious choice of PDE model being the semilinear problem 

u t = Lu + N(u), (17) 

where L is a linear differential operator and N(u) is some nonlinear map, typically con- 



taining derivatives of lower order than L. After discretizing in space, (17) is turned into a 
system of rid ODEs, for some large rid, L becomes an rid x n^-matrix, and TV" : R nd — > R nd a 
nonlinear function. We may now as in [52] introduce the action on R nd by some subgroup 
of the affine group represented as the semidirect product G = GL(rid) X R™ d . The group 
product, identity, and inverse are given as 

(A 1 ,b 1 )-{A 2 ,b 2 ) = {A 1 -A 2 ,A l b 2 + b l ), e = (J,0), (A, by 1 = (A~\ -A~ l b). 

The action on R™ d is 

(A, b)-x = Ax + b, (A, b)eG, xe R nd 
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and the Lie algebra and commutator are given as 

= (£,c), £€gl(n d ), cer, [(&, ci), (6, c 2 )] = ([6, 6], 6c 2 - 6ci + ci). 

In many interesting PDEs, the operator L is constant, so it makes sense to consider the 
rid + 1-dimensional subalgebra ql of g consisting of elements (aL, c) where a G R, c G R d , 
so that the map / : R n< * — >• is given as 

f(u) = (L,N(u)) 

One parameter subgroups are obtained through the exponential map as follows 

exp(t(L,6)) = (exp(tL),(j)(tL)tb). 

Here the entire function <j)(z) = (exp(z) — l)/z familiar from the theory of exponential 
integrators appears. As an example, one could now consider the Lie-Euler method in this 
setting, which coincides with the exponential Euler method 

u\ = exp(h(L, N(uq)) ■ uo = exp(hL)iiQ + hcf)(hL) N (uq) 

There is a large body of literature on exponential integrators, going almost half a century 
back in time, see |32j and the references therein for an extensive account. 



2.3 Isotropy — challenges and opportunities 

An issue which we have already mentioned a few times is that the map A* : — > X{M) 
defined in ([9]) is not necessarily injective. This means that the choice of / : M — > Q is not 
unique. In fact, if g: M — > q is any map satisfying \*(g(m))(m) = for all m £ M, then 



we could replace the map / by / + g in ( 10 ) without altering the vector field F. But such 
a modification of / will have an impact on the numerical schemes that we consider. This 
freedom in the setup of the methods makes it challenging to prove general results for Lie 
group methods, it might seem that some restrictions should apply to the isotropy choice 
for a more well defined class of schemes. However, the freedom can of course also be taken 
advantage of to obtain approximations of improved quality. 

An illustrative example is the two-sphere S 2 acted upon linearly by the special orthog- 
onal group SO(3). Representing elements of the Lie algebra so(3) by vectors in R 3 , and 
points on the sphere as unith length vectors in R 3 , we may facilitate (10) as 



F(m) = f(m) x m = (f(m) + a(m)m) x m, 

for any scalar function a: S 2 — > R. Using for instance the Lie-Euler method one would 
get 

mi = exp(/(m ) + a(m )m ) ■ m , (18) 
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where the exp is the matrix exponential of the 3x3 skew-symmetric matrix associated to 
a vector in R 3 via the "hat map". Clearly the approximation depends on the choice of 
a(m). The approach of Lewis and Olver |41j was to use the isotropy to improve certain 
qualitative features of the solution. In particular, they studied how the orbital error could 
be reduced by choosing the isotropy in a clever way. In Figure [2] we illustrate the issue of 
isotropy for the Euler free rigid body equations. The curve drawn from the initial point zq 
to Z\ is the exact solution, i.e. the momenta in body coordinates. The broken line shows 
the terminal points using the Lie-Euler method for a varying between and 25. 

Another potential difficulty with isotropy is the increased computational complexity 
when the group G has much higher dimension than the manifold M. This could for instance 
be the case with the Stiefel manifold St(d, k) if d ^> k. Linear algebra operations used in 
integrating differential equations on Stiefel manifold should preferably be of complexity 
0(dk 2 ). But solving a corresponding problem in the Lie algebra so(d) would typically 
require linear algebra operations of complexity 0(d 3 ). By taking advantage of the many 
degrees of freedom provided by the isotropy, it is actually possible to reduce the cost down 
to the required 0(dk 2 ) operations as explained for instance in [15] and [35] . 
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3 Symplectic integrators on the cotangent bundle of a Lie 
group 

In this section we shall assume that the manifold is the cotangent bundle T*G of a Lie 
group G. Let R g : G — > G be the right multiplication operator such that R g {h) = hg for any 
h G G. The tangent map of R g is denoted R g * := TR g . Any cotangent vector p g G T*G 
can be associated to [i G g* by right trivialization as follows: Write v g G T g G in the form 
v g = R g *£, where £ G g, so that (p g ,v g ) = (p g , R g *£,) = (R*p g ,£,), where we have used R* for 
the dual map of R g *. We therefore represent p g G T*G by fi = R* g p g G 0*. Thus, we may 
use as phase space G x g* rather than T*G. For applying Lie group integrators we need a 
transitive group action on G x g* and this can be achieved by lifting the group structure 
of G and use left multiplication in the extended group. The semi-direct product structure 
on G := G k g* is defined as 

{gi, Mi) • {92, 1^2) = {gi ■ g2, \n + Ad*_i^ 2 ) (19) 

Similarly, the tangent map of right multiplication extends as 

TR(g, fl )(R h *C,^) = (Rhg* 0-ad£Ad£_i/x), g,heG, ( G g, /u,^Gg*. 
Of particular interest is the restriction of TR/ g ^ to T e G = g k g*. 

The natural symplectic form on T*G is defined as 

^( g ,p g )((M,<57ri), (8v2,6ir 2 )) = {5tt 2 ,5vi) - (5iri,5v 2 ), 
and by right trivialization it may be pulled back to G ix g* and then takes the form 

u >(g,n)(( R g*£i> Sv i)>( R g*£2> Sv 2)) = (<^2,£i) - (<^i,6) - (M6,6]> (20) 

The presentation of differential equations on T*G is now achieved via the action by left 
multiplication, meaning that any vector field F G X(G ix g*) is expressed by means of a 
map /: G k g* g k g* 

F( 5 ,/i) = T e R M f(g,fi) = (R g *h,f 2 - ad^/x) (21) 

where f\ = fi(g,/J.) G g, f 2 = f2(g,l^) G g* are the two components of /. We are 
particularly interested in the case that F is a Hamiltonian vector field which means that 
F satisfies the relation 

\ f oj = dH (22) 
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for some Hamiltonian function H : T*G — > R. A simple calculation using (20), (21) and 



(22) shows that the corresponding map / for such a Hamiltonian vector field is 



We have come up with the following family of symplectic Lie group integrators on G x g* 
(&,rii) = hf(Gi,Mi), n* = Ad* xp(x .)nj, i = l,...s. 



(9i,Vi) = exp(T, dexp/* ^ 6^) • (#o,Mo) 



i=l 



y = ^6^i, ,Y, ^",.,s ; . i = l,...,s. 

t=l j=l 

Gi = exp(Xi) ■ g , i = l,...,s 

s ha- 
Mi = dexply(/i) + > (6jdexply ^-^dexp^x )nj , i = 1, . . . , s. 

3=1 bl 

It is assumed that Ylt=i h = I and that 6, / 0, 1 < i < s. The symplecticity of these 
schemes is a consequence of their derivation from a variational principle, following ideas 
similar to that of [3] and [61]. One should be aware that order barriers for this type of 
schemes may apply, and that further stage corrections may be necessary to obtain high 
order methods. In Figure [3] we show numerical experiments for the heavy top where the 
Hamiltonian is given as 

where I: q — > q* is the inertia tensor, here represented as a diagonal 3x3 matrix, uq is 
the initial position of the top's center of mass, and e3 is the canonical unit vector in the 
vertical direction. We have chosen I = 10 3 diag(l, 5, 6) and uq = e3. The initial values used 
were go = I (identity matrix) and no = 101 (1, 1, 1) T . We have used the coefficients of the 
well known ^-method, where s = 1, b\ = 1, a\\ = 8, and we compare the behaviour of the 
symplectic schemes presented here to the Runge-Kutta-Munthe-Kaas (RKMK) method 
with the same 6- value. In Figure [3] we have drawn the time evolution of the center of 
mass, u n = g n • uo, the characteristic band structure observed for the symplectic methods 
was reported in |17| . The RKMK method with 6 = ^ exhibits a similar behaviour, but 
the bands are expanding faster than for the symplectic ones. We have also found in these 
experiments that none of the symplectic schemes, 6 = and 6 = \ have energy drift, but 
this is also the case for the RKMK method with = ^. For 9 = 0, however, the RKMK 
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method shows energy drift as expected. These tests were done with step size h = 0.05 over 
10 5 steps. 



4 Discrete gradients and integral preserving methods on Lie 
groups 

The discrete gradient method for preserving first integrals has to a large extent been made 
popular through the works of Gonzalez [27] and McLachlan et al. [46] . the latter proved 
the result that under relatively general circumstances, a differential equation which has a 
first integral I(x) can be written in the form 

x = S{x)VI(x) 

for some non- unique solution-dependent skew-symmetric matrix S(x). The idea is to intro- 
duce a mapping which resembles the true gradient; a discrete gradient VI: R d x R d — > H d 
is a continuous map which satisfies the following two conditions 

Vl(x,x) = W(x), Vx, 

I(y)-I(x) = V I(x,y)(y-x), Vx^y. 
An integrator which preserves /, that is, I(x n ) = I(xq) for all n is now easily devised as 



h 



S(x n , x n+1 )VI(x n , x n+ i) 



where S(x,y) is some consistent approximation to S(x), i.e. S(x,x) = S(x). There exist 
several discrete gradients, two of the most popular are 

VJ(z, y) = I VI((y + (1 - ()x) d( (23) 
J o 



and 



VI(x,y) = VI[— + ||y-g|| 2 (y~ x )- ( 24 ) 



The matrix S(x, y) can be constructed with the purpose of increasing the accuracy of the 
resulting approximation, see e.g. [59] . 

We now generalize the concept of the discrete gradient to a Lie group G. We consider 
differential equations which can, for a given dual two-form^] to G 02(G) and a function 
H : G — > R be written in the form 

x = i dH u (25) 



x By dual two- form, we here mean a differential two- form on G such that on each fiber of the cotangent 
bundle we have to\ x : T*G x T*G — > R, a bilinear, skew-symmetric form. 
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where i a is the interior product i a u((3) = u(a, f3) for any two one-forms a, (3 £ J7 1 (G ! ). The 
function H is a first integral since 

j t H{x(t)) = dH\ x{t) {x(t)) = u(dH,dH) = 0. 

We define the trivialized discrete differential (TDD) of the function H to be a continuous 
map dH : G x G — >■ Q* such that 

H(x')-H(x) = $H{x,x'),\og{x'x- 1 )) 

d~H(x,x) = R* x dH\ x 

A numerical method can now be defined in terms of the discrete differential as 

x = exp(/i i^ H ^ x y)^(x, x')) x 

where u) : G x G — > A 2 (g*) is continuous. This exterior form is some local trivialized 
approximation to u. For consistency we require 

u)(x,x)(R* x a,R* x P) = uj\ x (a,f3), a,(3£T*G. (26) 

We easily see that this method preserves H exactly, since 

H{x) - H(x) = (dtf^x'yog^'ar 1 )) 

= (dH{x,x'),hi^ H{xxl) u;) 

= ho(dH(x, x'),dH(x, x')) = 0. 



Extending ( 23 ) to the Lie group setting we define the following TDD 

dH(x, x') = J R* m dH\ m dZ, £(0 = exp(£ log(x' • x'^x. (27) 
Similarly, for any given inner product on q, we may extend the discrete gradient (24) to 
dH(x,x') = R%dH\, + H{X ' ] ~ H{X } - {R%m ^\ \ V = log(x'x-i), (28) 

WvW 

where i£G could for instance be x = exp(r//2)x, a choice which would cause dH(x, x') = 
d~H{x',x). 
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Suppose that the ODE vector field F is known as well as the invariant H. A dual 
two-form uj can now be defined in terms of a Riemannian metric on G. By index raising 
applied to dH, we obtain the Riemannian gradient vector field gradif, and we define 



grad H A F 
||grad#|| 2 



i dH uj = F. 



Example: We consider the equations for the attitude rotation of a free rigid body- 
expressed using Euler parameters. The set S* 3 = {q £ R 4 | ||<q||2 = 1} with q = [go>q] 
(qo G R and q € R 3 ), is a Lie group with the quaternion product 

p • q = \poqo - p T q, pop + <?oq + p x q] 

with unit e = [10 0] and inverse q c = [qo, — q]. Denote by " * " the hat-map identifying 
R 3 with so (3): 





VI 







—V3 


V2 


V = 


V2 


!->■ V = 


^3 





-Vl 




V3 




— V2 


VI 






The Lie group can be mapped into 50(3) by the Euler-Rodrigues map: 

£(q) =/ 3 + 2g q + 2q 2 , 

where /3 denotes the 3x3 identity matrix. The Lie algebra s 3 of S* 3 is the set of so called 
pure quaternions, the elements of R 4 with first component equal to zero, identifiable with 
R 3 and with so (3) via the hat-map. 

The equations or the attitude rotation of a free rigid body on S 3 read 



and 



q = /(q) • q, /(q) = q • v • q c , 



v = [0 , v], v = - I 1 <5(q c )m , 



mo is the initial body angular momentum and I is the diagonal inertia tensor, and according 
with the notation previously used in this section F(q) = /(q) • q. The energy function is 

#(q) = - m„ £:(q)I _1 <S(q c )m . 

We consider the R 3 Euclidean inner product as metric in the Lie algebras 3 , and obtain 
by right translation a Riemannian metric on S 13 . The Riemannian gradient of H with 
respect to this metric is then 



grad#= (h ~ qq T )V#, 
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where I4 is the identity in R 4 and VH is the usual gradient of H as a function from R 4 
to R. We identify s 3 with its dual, and using gradH in (24) we obtain the (dual) discrete 



differential grad-£f(q, q') 6 s 3 . 

The two-form uj = ^gradJJ|p w ^h respect to the right trivialization can be identified 
with the 4x4 skew-symmetric matrix 



Wi?(q) = jj— , 7 £ « . ^ • q = -^(q), 7 • «1 = gradi^(q), 



cj^(q) has first row and first column equal to zero. We choose oj to be 

w(q, q') = Wfl(q), q = exp(r?/2)q, 77 = log(q • q c ), 

i.e. lur frosen at the mid-point q. The energy-preserving Lie group method of second order 
is 



q' = exp(/icj(q, q') grad H (q, q')) • q, 

and exp is the exponential map from s 3 to S 3 , log: S s — > s 3 is its inverse. 

In Figure[4]we plot the body angular momentum vector m = £ (q c )mo on a time interval 
[0, T], T = 1000, for four different methods: the Lie group energy-preserving integrator 
just described (top left), the built-in Matlab routine ode45 with absolute and relative 
tolerance 10 -6 (top right); the ode45 routine with tolerances 10 -14 (bottom left); and the 
explicit Heun RKMK Lie group method (bottom right). The two Lie group methods both 
have order 2. The energy preserving method is both symmetric, energy preserving and it 
preserves the constraint ||q||2 = 1- The Lie group integrators use a step-size h = 1/64. 
The solution of the built-in Matlab routine at high precision is qualitatively similar to 
the highly accurate solution produced by Matlab with tolerances 10 -14 , the energy error 
is also comparable for these two experiments. The performance of other Matlab built-in 
routines we tried was worse than for ode45. We remark that the equations are formulated as 
differential equations on S 3 , a formulation of the problem in form of a differential algebraic 
equation would possibly have improved the performance of the Matlab built-in routines. 
However it seems that the preservation of the constraint alone can not guarantee the good 
performance of the method. In fact the explicit (non-symmetric) Lie group integrator 
preserves the constraint ||q||2 = 1, but performs poorly on this problem (see Figure [4] 
bottom right). The cost per step of the explicit Lie group integrator is much lower than 
for the energy-preserving symmetric Lie group integrator. 



5 Applications to nonlinear problems of evolution in classi- 
cal mechanics 

The emphasis on the use of Lie groups in modeling and simulation of engineering problems 
in classical mechanics started in the eighties with the pioneering and fundamental work of 
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Figure 4: Free rigid body angular momentum, time interval [0, 1000], moments of inertia 
I\ = 1, I2 = 5, I3 = 60, initial angular velocity I mo = [1, 1/2, —1]: (top left) energy- 
preserving Lie group method, h = 1/64; (top right) ode45 with tolerances 10 -6 ; (bottom 
left) ode45 with tolerances 10~ 14 ; (bottom right) Heun RKMK, h = 1/64. 
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Figure 5: Geometric rod model. 92 line of centroids, cross section identified by the frame 
[ti,t 2 ,t 3 ]. 

J.C. Simo and his collaborators. In the case of rod dynamics, for example, models based 
on partial differential equations were considered where the configuration of the centreline 
of the rod is parametrized via arc-length, and the movement of a rigid frame attached to 
each of the cross sections of the rod is considered, see Figure [5] This was first presented 
in a geometric context in [64J. 

In robot technology, especially robot locomotion and robot grasping, the occurrence of 
non-holonomically constrained models is very common. The motion of robots equipped 
with wheels is not always locally controllable, but is often globally controllable. A classical 
example is the parking of a car that cannot be moved in the direction orthogonal to its 
wheels. The introduction of Lie groups and Lie brackets to describe the dynamics of such 
systems, has been considered by various authors, see for example |56j. The design of 
numerical integration methods in this context has been addressed in the paper of Crouch 
and Grossman, |20j. These methods have had a fundamental impact to the successive 
developments in the field of Lie group methods. 

The need for improved understanding of non-holonomic numerical integration has been 
for example advocated in |47j . Recent work in this field has led to the construction of 
low order non-holonomic integrators based on a discrete Lagrange-d'Alembert principle, 
[T9"l 14*9"] . The use of Lie group integrators in this context has been considered in [101 00] . 

We have already mentioned the relevance of rigid body dynamics to the numerical 
discretization of rod models. There are many other research areas in which the accurate 
and efficient simulation of rigid body dynamics is crucial: molecular dynamics, satellite 
dynamics, and celestial mechanics just to name a few, [39]. In some of these applications, 
it is desirable to produce numerical approximations which are accurate possibly to the 
size of roundoff. The simulations of interest occur over very long times and/or a large 
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number of bodies and this inevitably causes propagation of errors even when the integrator 
is designed to be very accurate. For this reason accurate symplectic rigid body integrators 
are of interest because they can guarantee that the roundoff error produced by the accurate 
computations can stay bounded also in long time integration. This fact seems to be of 
crucial importance in celestial mechanics, [37] • A symplectic and energy preserving Lie 
group integrator for the free rigid body motion was proposed in [12] • The method computes 
a time re-parametrization of the exact solution. Some recent and promising work in this 
field has been presented in |48[ [T71 [TU1 [29] . The control of rigid bodies with variational Lie 
group integrators was considered in [40J. 

In the next section we illustrate the use of Lie group methods in applications on a 
particular case study, the pipe-laying process from ships to the bottom of the sea. 

5.1 Rigid body and rod dynamics 

5.1.1 Pipe-laying problem 

The simulation of deep-water risers, pipelines and drill rigs requires the use of models of 
long and thin beams subject to large external forces. These are complicated nonlinear 
systems with highly oscillatory components. We are particularly interested in the correct 
and accurate simulation of the pipe-laying process from ships on the bottom of the sea, 
see Figure [6j The problem comprises the modeling of two interacting structures: a long 
and thin pipe (modeled as a rod) and a vessel (modeled as a rigid body). The system is 
subject to environmental forces (such as sea and wind effects). The control parameters for 
this problem are the vessel position and velocity, the pay-out speed and the pipe tension 
while the control objectives consist in determining the touchdown position of the pipe as 
well as ensuring the integrity of the pipe and to avoid critical deformations, [2H [62] . 

The vessel rigid body equations are determining the boundary conditions of the rod. 
They can be integrated numerically with a splitting and composition technique where the 
vessel equations are split into a free rigid body part and a damping and control part. The 
free rigid body equations can be solved with a method proposed in [T7] where the angular 
momentum is accurately and efficiently computed by using Jacobi elliptic functions, the 
attitude rotation is obtained using a Lie group method based on the Magnus expansion, 
and the control and damping part are solved exactly. 

Simulations of the whole pipe-lay problem with local parametrizations of the pipe and 
the vessel based on Euler angles have been obtained in [53] . 

5.1.2 Rod dynamics 

At fixed time each cross section of the pipe is the result of a rigid rotation in space of 
a reference cross section, and, analogously, for each fixed value of the space variable the 
corresponding cross section evolves in time as a forced rigid body, see Figure [5] 
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Pipeline vessel 




Figure 6: The pipe-lay process. 



In [66J, partitioned Newmark integrators, of Lie group type, and of moderate order were 
considered for this problem. While classical Newmark methods are variational integrators 
and as such are symplectic when applied to Hamiltonian systems |45] . the particular ap- 
proach of [66j involves the use of exponentials for the parametrization of the Lie group 
SO(3) and, the symplecticity in this case is not easily achieved. Moreover since the model 
is a partial differential equation space and time discretizations should be designed so that 
the overall discrete equations admit solutions and are stable. It turns out that conven- 
tional methods perform poorly on such problem in long time simulations. To obtain stable 
methods reliable in long-time simulation, an energy-momentum method was proposed for 
the rod problem in [65J. Later this line of thought has been further developed in |60j. 
The Hamiltonian formulation of this model allows to derive natural structure preserving 
discretizations into systems of coupled rigid bodies |43j . 

Following the geometric space-time integration procedure proposed in [26], a multi- 
Hamiltonian formulation^] of these equations has been proposed in [18], using the Lie group 
of Euler parameters. The design of corresponding multi-symplectic Lie group discretiza- 
tions is still under investigation. 

2 For a definition of the multi-symplectic structure of Hamiltonian partial differential equations, see [4]. 
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6 Applications to problems of data analysis and statistical 
signal processing 

The solution of the many-body Schrodinger eigenvalue problem, 

HI' = EV, (29) 

where the so called electronic ground state (the smallest eigenstate) is sought, is an impor- 
tant problem of computational chemistry. The main difficulty is the curse of dimensionality. 



Since H is a differential operator in several space dimensions, a realistic simulation of (29) 
would require the numerical discretization and solution of a partial differential equation 
in several space dimensions. The number of space dimensions grows with the number of 
electrons included in the simulation. 

The eigenvalue problem admits an alternative variational formulation. Instead of look- 
ing for the smallest eigenvalue and eigenfunction of the Schrodinger equation, one minimizes 
directly the ground state energy. 

After appropriate spatial discretization, the problem becomes a minimization problem 
on a Riemannian manifold M, 

min (/>(#), (30) 
xeM 

where <j) : M. — > R is a smooth discrete energy function to be minimized on M. The discrete 
energy considered here is the so called Kohn-Sham energy, pQ . For a related application 
of Lie group techniques in quantum control see |21| . 



The general optimization problem giving rise to (30) arises in several applied fields 
ranging from engineering to applied physics and medicine. Some specific examples are 
principal component/subspace analysis, eigenvalue and generalized eigenvalue problems, 
optimal linear compression, noise reduction, signal representation and blind source sepa- 
ration. 

6.1 Gradient-based optimization on Riemannian manifolds 

Let Ai be a Riemannian manifold with metric (-,-)m an d 4>- M. — > R be a smooth cost 



function to be minimized on A4, and we want to solve (30). The optimization method 
based on gradient flow, (written for the minimization problem only, for the sake of easy 
reading), consists in setting up the differential equation on the manifold, 

x(t) = -grad(p(x(t)), (31) 



with appropriate initial condition x(0) = xq £ JM. The equilibria of equation (31) are 
the critical points of the function <p. In the above equation, the symbol grad</> denotes 
the Riemannian gradient of the function <fi with respect to the chosen metric. Namely, 
grad(p(x) G T X M and <f>'\ (v) = (grad^(x), v)m f° r an v £ T X M. 
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The solution of (31 ) on Ai may be locally expressed in terms of a curve on the tangent 
space T Xo A4 using a retraction map (p. Retractions are tangent space parametrizations of 
Ai, and allow us to write 

x(t) = <p X0 (a(t)), a(t) e T X0 M, t G [0,*/], 

for small enough tf, see [63] for a precise definition. 

In most applications of interest, see for example [5j[3T], Ai is a matrix manifold endowed 
with a Lie group action and there is a natural way to define a metric and a retraction. In 
fact, let At be a manifold acted upon by a Lie group G, with a locally transitive group 
action A(g,x) = A. x (g). Let us also consider a coordinate map tp, 

ij>: 0^G, and a, := (A x o^)'\ Q . (32) 

One can prove that if there exists a linear map & x : T x Ai — > g such that p x o a x = Id^ m j 
then cp x , given by 

^(v) -(Aj-o^oaj;)^), (33) 

is a retraction, see [16]. The existence of a x is guranteed, at least locally, by the transitivity 
of the action and the fact that ip is a local diffeomorphism. The approach is analogous 
to the one described for differential equations in section [2~T| Therefore, we can construct 
retractions using any coordinate map from the Lie algebra, g, to the group. 
Any metric on q, (•, -) induces a metric on Ai by 

(v x ,w x ) M = (a>x(v x ),a x (w x )) s . 

Also, we may define the image of the tangent space under the map & x : 

m x := & X (T X M) C 0. (34) 

The set is a linear subspace of the Lie algebra q, often of lower dimension. Parametriza- 
tions of the solution of (31) involving the whole Lie algebra are in general more computa- 
tionally intensive than those restricted to m^, but, if the isotropy is chosen suitably, they 
might lead to methods which converge faster to the optimum. 

For the sake of illustration, we consider the minimization on a two-dimensional torus 
T 2 = S 1 x S . Here we denote by S 1 the circle, i.e. 

S 1 = {g-e 1 \ R 2 , g(a)GSO(2)}, 



g(a) = exp(aE), E 



-1 

1 



< a < 2vr, 



where ei is the first canonical vector and 50(2) is the commutative Lie group of planar 
rotations. Any element in T 2 is of the form 

x G T 2 , x = (g(e)e 1 ,g(i P )e 1 ), g(0), g(<p) G 50(2). 
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The Lie group acting on T 2 is 50(2) x 50(2), its corresponding Lie algebra is so(2) xso(2), 
which has dimension d = 2 and basis {(E, O), (O, E)}, where O is the zero element in so(2). 
The Lie group action is 

AvofaM = (hi ■ g(9)ei, h x ■ g(<p)ei), (h u h 2 ) G 50(2) x 50(2), 

and ip = exp. Any v Xo G T XQ T 2 can be written as 

v Xo = (aEe 1 ,/3Eei) 

for some a, f3 G R, so 

a>x (vxo) = (aE,(3E). 

Assume the cost function we want to minimize is simply the distance from a fixed plane 
in R 3 , say the plane with equation y = 8. This gives 

<j,{g(e)e 1 ,g( t p)e 1 ) = |(1 + cos(0)) sin(^) - 8|, (35) 

and the minimum is attained in 9 = and <~p = In Figure [7] we plot — grad0, the 

negative gradient vector field for the given cost function. The Riemannian metric we used 
is, 

((aiEe!, PiEei), (a 2 Eei, /3 2 Eex))) T 2 = a x a 2 + 

and {a.\Ee\, (3\Eei) G T( ei ei )T 2 . This metric can be easily interpreted as a metric on the 
Lie algebra q = so (2) x so (2): 

((axE, faE), (a 2 E, f3 2 E)) Q = a x a 2 + Pifo- 

At the point po = (g(6o)ei, g(tpo)ei) G T 2 = 5 1 x 5 1 the gradient vector field can be 
represented by 

(7 • E ■ g(6 )ei,S ■ E • gfafiei), 
where 7 and 5 are real values given by 

7 = -C ■ sin(6>o) sm((p ), d = C • (1 + 003(6*0)) cos(</? ), 

and 

C = 2((1 + cos(0 o ) sin(v? ) - 8). 

Gradient flows are not the only type of differential equations which can be used to 
solve optimization problems on manifolds. Alternative equations have been proposed in 
the context of neural networks [T2] . Often they arise naturally as the Euler-Lagrange 
equations of a variational problem. 

3 We have used a parameterization of T 2 in R 3 in angular coordinates, obtained applying the following 
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Figure 7: The gradient vector field of the cost function (f>((g(9)ei, g((p)e±) = (1 + eJ g(9)ei)- 
e^g((^)ei — 8) 2 on the torus. The vector field points towards the two minima, the global 
minimum is marked with a black spot in the middle of the picture. 
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6.2 Principal component analysis 



Data reduction techniques are statistical signal processing methods that aim at providing 
efficient representations of data. A well-known data compression technique consists of 
mapping a high-dimensional data space into a lower dimensional representation space by 
means of a linear transformation. It requires the computation of the data covariance 
matrix and then the application of a numerical procedure to extract its eigenvalues and 
the corresponding eigenvectors. Compression is then obtained by representing the signal in 
a basis consisting only of those eigenvectors associated with the most significant eigenvalues. 

In particular, Principal Component Analysis (PCA) is a second-order adaptive statis- 
tical data processing technique that helps removing the second-order correlation among 
given random signals. Let us consider a stationary multivariate random process x(t) G R n 
and suppose its covariance matrix A = E[(x — E[x])(x — E[x])] T ] exists and is bounded. 
Here the symbol E[-] denotes statistical expectation. If A G R nx ™ is not diagonal, then 
the components of x(t) are statistically correlated. One can remove this redundancy by 
partially diagonalizing A, i.e. computing the operator F formed by the eigenvectors of the 
matrix A corresponding to its largest eigenvalues. This is possible since the covariance ma- 
trix A is symmetric (semi) positive-definite, and F G St(n,p). As a consequence the new 
random signal defined by y{t) := F T (x(t) — E[x(t)]) G R p has uncorrelated components, 
with p < n properly selected. The component signals of y(t) are the so called principal 
components of the signal x(t), and their relevance is proportional to the corresponding 
eigenvalues af = E[yf] which here are arranged in descending order (<r| > o~f + 1 ). 

Thus, the data stream y(t) is a compressed version of the data stream x(t). After the 
reduced-size data has been processed (i.e. stored, transmitted), it needs to be recovered, 
that is, it needs to be brought back to the original structure. However, the principal- 
component-based data reduction technique is not lossless, thus only an approximation 
x(t) G R n of the original data stream may be recovered. As F is a tall-skinny orthogonal 
operator, an approximation of x{t) is given by x(t) = Fy(t) + E[x\: Such approximate data 
stream minimizes the reconstruction error E\\\x — which equals ^?=n+i a i' 

For a scalar or a vector-valued random variable x G R n endowed with a probability 
density function p x : x G R n — > p x (x) G R, the expectation of a function (3: R n — > R is 
defined as 



Under the hypothesis that the signals whose expectation is to be computed are ergodic, the 



with < 6, (p < 2-k. This is equivalent to the composition of two planar rotations and one translation in 




(36) 



mapping 




R 3 . 
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actual expectation (ensemble average) may be replaced by temporal-average on the basis 
of the available signal samples, namely 

£[fl«if>(s(t)). (37) 
6.3 Independent Component Analysis 

An interesting example of a problem that can be tackled via statistical signal processing is 
the cocktail-party problem. Let us suppose n signals xi(t), . . . ,x n (t) were recorded from n 
different positions in a room where there are p sources or speakers. Each recorded signal 
is a linear mixture of the voices of the sources si(t), . . . , s p (t), namely 

xi(t) = ai^si (£) H h a 1>p s p (t), 



x n (t) = a„,isi(i) H h a ntP s p (t), 

where the n-p coefficients aij £ R denote the mixing proportions. The mixing matrix A = 
(cLij) is unknown. The cocktail party problem consists in estimating signals s\(t),. . . , s p (t) 
from the only knowledge of their mixtures xi(t), . . . , x n (t). The main assumption on the 
source signals is that s±(t) , . . . , s p (t) are statistically independent. This problem can be 
solved using Independent Component Analysis. 

Typically, one has n > p, namely, the number of observations exceeds the number of 
actual sources. Also, a typical assumption is that the source signals are spatially white, 
which means E[ss T ] = I p . The aim of independent component analysis is to find estimates 
y(t) of signals in s(t) by constructing a de-mixing matrix W G R nx P an d by computing 
y(t) := W T x(t). Using statistical signal processing methods, the problem is reformulated 
into an optimization problem on homogeneous manifolds for finding the de-mixing matrix 
W. 

The geometrical structure of the parameter space in ICA comes from a signal pre- 
processing step named signal whitening, which is operated on the observable signal x(t) — > 
x(t) £ R p in such a way that the components of the signal x(t) are uncorrelated and have 
variances equal to 1, namely E[xx ] = I p . This also means that redundant observations 
are eliminated and the ICA problem is brought back to the smallest dimension p. This 
can be done by computing E[xx T ] = VDV T , with V G St(n,p) and D G R pxp diagonal 
invertible. Then pre-whitening is obtained as 

x(t) := D~W T x{t), (38) 
and E[xx T ] = AE[ss T ]A T = AA T = I p . 
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After observable signal pre-whitening, the de-mixing matrix may be searched for such 
that it solves the optimization problem 



max 6{W). (39) 

As explained, after pre-whitening, the number of projected observations in the signal 
x(t) equals the number of sources. However, in some applications it is known that not all 
the source signals are useful, so it is sensible to analyze only a few of them. In these cases, 
if we denote by p <C p the actual number of independent components that are sought after, 
the appropriate way to cast the optimization problem for ICA is 

max 4>(W), with p <C p. (40) 

WGSt(n,p) 

As a possible principle for reconstruction, the maximization or minimization of non- 
Gaussianity is viable. It is based on the notion that the sum of independent random 
variables has distribution closer to Gaussian than the distributions of the original random 
variables. A measure of non-Gaussianity is the kurtosis, defined for a scalar signal z G R 
as 

kurt(z) := E[z*] - 3E 2 [z 2 ]. (41) 

If the random signal z has unitary variance, then the kurtosis computes as kurt(z) = E[z 4 ] — 
3. Maximizing or minimizing kurtosis is thus a possible way of estimating independent 
components from their linear mixtures, see [12] and references therein for more details. 



6.4 Computation of Lyapunov exponents 

The Lyapunov exponents of a continuous dynamical system x' = F(x), (x(t) £ M. n ) provide 
a qualitative measure of its complexity. They are numbers related to the linearization A(t) 
of x' = F{x) along a trajectory x(t). Consider the solution U of the matrix problem 

U = A(t)U, 17(0) = U , n x n. 

The logarithms of the eigenvalues of the matrix 

A = lim (U(t) T U(t))^ , 

are the Lyapunov exponents for the given dynamical system. In [22J a procedure for 
computing just k of the n Lyapunov exponents of a dynamical system is presented. The 
exponents are computed by solving an initial value problem on St(n, k) and computing 
a quadrature of the diagonal entries of a k x k matrix-valued function. The initial value 
problem is defined as follows: 

Q= (A- QQ T A + QSQ T ) Q, 
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with random initial value in St(n, k) and 



Sk, 



(Q 1 AQ) kJ , 
0, 

-(Q T AQ)j k , 



k>j, 
k = j, 
k<j, 



k,j = l,...,p. 



It can be shown that the i-th. Lyapunov exponent Aj can be obtained as 



lim - 

t— >co t 



Bii(s)ds, 



1, . . . , k, 



(42) 



and 



B = Q AQ - S. 



One could use for example the trapezoidal rule to approximate the integral ( 42 ) and com- 
pute Aj (i = 1, . . . , k). We refer to [22] for further details on the method, and to |14j for 
the use of Lie group integrators on this problem. Lie group methods for ODEs on Stiefel 
manifolds have also been considered in [THl ESI H3] • 



We have here presented a selection of applications of Lie group integrators with par- 
ticular emphasis to problems of evolution in classical mechanics and problems of signal 
processing. This is by no means an exhaustive survey; other interesting areas of applica- 
tion are for example problems in vision and medical imaging. 
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